## b (i.e., beta) is the null hypothesis to test.
#t.r - instrumented outcome
#t.d - 
#c.r - nonintrumented outcome
#c.d -

# Luke uses now
# iv_pval <- function(t.r, t.d, c.r, c.d, b) {
#   z <- (t.r - b*t.d) - (c.r - b*c.d)
#   n <- length(z)
#   r <- rank(abs(z))
#   tplus <- sum(r[z > b])
#   lower <- psignrank(tplus, n, lower.tail=FALSE)
#   upper <- psignrank(tplus, n, lower.tail=TRUE)
#   min(lower,upper)
# }

HodgesLehman <-
function(Rt, Rc, Dt, Dc, Zt = rep(1, length(Rt)),
         Zc = rep(0, length(Rt)), b = 0, alpha = 0.025,
         Beta = 2, BetaInc = 0.01) {
  
  ## Length of the response vector for the treated.
  n <- length(Rt)
  
  ## Check that the response and dose vectors have the same length.
  if (length(Rc) != n | length(Dt) != n | length(Dc) != n |
    length(Zt) != n | length(Zc) != n) {
    stop("Length of response, dose, and encouragement vectors not equal.")
  }

  iv_pval <- function(b = 0, return_tplus = FALSE) {
    z <- (Rt - b*Dt) - (Rc - b*Dc)
#     z <- z[Rt - Rc != 0]
    n <- length(z)
    r <- rank(abs(z))
    tplus <- sum(r[z > b])
    if(return_tplus) return(tplus)
    lower <- psignrank(tplus, n, lower.tail = FALSE)
    upper <- psignrank(tplus, n, lower.tail = TRUE)
    return(min(lower, upper))
  }
  
  pval <- iv_pval(b = 0, return_tplus = FALSE)
  
  # Confidence interval
  rng <- seq(-abs(Beta), abs(Beta), by = BetaInc)
  ci <- sapply(rng, FUN = function(x) iv_pval(b = x, return_tplus = FALSE))
  ci.out <- data.frame(Beta = rng, pval = round(ci, 3))
  
  ## Estimate the p-value for all null hypothesis in range as
  ## defined by Beta incremented by BetaInc.
#   rng <- seq(-abs(Beta), abs(Beta), by = BetaInc)
#   sens <- rep(NA, length(rng))
#   for (i in 1:length(rng))
#     sens[i] <- iv_pval(Rt, Rc, Dt, Dc, b = rng[i])
#   Sens <- data.frame("Beta" = rng, "pval" = sens)
  
  ## Select some reasonable intervals from which to begin the
  ## search. Supress warnings, which will be dealt with cleanly below.
  x2 <- min(ci.out$Beta[ci.out$pval > alpha])
  x3 <- x2 + 2*BetaInc; x1 <- x2 - 2*BetaInc

  y2 <- max(ci.out$Beta[ci.out$pval > alpha])
  y3 <- y2 + 2*BetaInc; y1 <- y2 - 2*BetaInc

  ## This is, in effect, a cost function that will be passed to
  ## optimize() to minimize.
  maim_sens <- function(b) {
    abs(iv_pval(b) - alpha)
  }
  
  ## Find the lower and upper bounds for beta only if boundaries were
  ## found from the initial sensitivity estimates.
  if (x2 == -Beta | any(is.infinite(c(x1, x2, x3)))) {
    cat("Lower boundary undefined. Try increasing Beta.\n")
    b1 <- -Inf
  }
  else {
#     b1 <- iv_gssearch(x1, x2, x3, maim_dist, Rt, Rc, Dt, Dc, Zt, Zc)
#       b1 <- .gssearch(iv_pval, x1, x2, x3)
      b1 <- optimize(maim_sens, lower = x1, upper = x3)$minimum
                     
  }
  
  if (y2 == Beta | any(is.infinite(c(y1, y2, y3)))) {
    cat("Upper boundary undefined. Try increasing Beta.\n")
    b2 <- Inf
  }
  else {
#     b2 <- iv_gssearch(y1, y2, y3, maim_dist, Rt, Rc, Dt, Dc, Zt, Zc)
    b2 <- optimize(maim_sens, lower = y1, upper = y3)$minimum
  }
  
  ConfInt <- data.frame(lower = b1, upper = b2)
  rownames(ConfInt) <- c("Interval")

  ## Calculate Hodges-Lehmann point estimate. This should be the Beta
  ## for which the p-value is maximized. Where more than one point is
  ## returned, the mean midpoint is chosen (see Rosenbaum 2010,
  ## p. 136). If one of the bounds is not defined, a warning is
  ## returned and the HL estimate is set to Inf.
  if (any(is.infinite(c(b1, b2)))) {
    cat("Confidence interval not defined. Hodges-Lehmann estimate not calculated.\n")
    HL <- Inf
  }
  else {
    I <- length(Zt)
    T.null <- (I*(I+1))/4
    tplus <- sapply(rng, FUN = function(x) iv_pval(b = x, return_tplus = TRUE))
    HL <- rng[which.min(abs(tplus - T.null))]
  }
  
  Obj <- list("HL" = HL, "Alpha" = alpha, "Beta" = Beta,
              "BetaInc" = BetaInc, pval = pval,
              "Interval" = ConfInt)
  
  class(Obj) <- c("ivsens", class(Obj))
  Obj
}
